{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rate-control Empirical Analysis"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "<script>\n",
    "  function code_toggle() {\n",
    "    if (code_shown){\n",
    "      $('div.input').hide('500');\n",
    "      $('#toggleButton').val('Show Code')\n",
    "    } else {\n",
    "      $('div.input').show('500');\n",
    "      $('#toggleButton').val('Hide Code')\n",
    "    }\n",
    "    code_shown = !code_shown\n",
    "  }\n",
    "\n",
    "  $( document ).ready(function(){\n",
    "    code_shown=false;\n",
    "    $('div.input').hide()\n",
    "  });\n",
    "</script>\n",
    "<form action=\"javascript:code_toggle()\"><input type=\"submit\" id=\"toggleButton\" value=\"Show Code\"></form>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple linear regression\n",
    "\n",
    "We performed a simple linear regression of the bitrate with respect to the quantizer,\n",
    "operating on the logarithm of both.\n",
    "The data set used was all of the video clips on https://media.xiph.org/video/derf/\n",
    "as well as subset3 (for extra I-frame data).\n",
    "To enable processing an arbitrarily large data set, an online regression algorithm was implemented.\n",
    "In practice, [440MB of text formatted data](https://ba.rr-dav.id.au/data/rav1e/rc-data.tar.xz) were sufficient.\n",
    "\n",
    "The raw final state of the online regression for each segment follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{0: (2.7695336845023429016e+17,\n",
      "     -1.07133222745900214305e+17,\n",
      "     8.397389236899963756e+38,\n",
      "     1.1092601296011764081e+39,\n",
      "     91218,\n",
      "     9581914802246888,\n",
      "     419495824559273192),\n",
      " 1: (2.5070212498929263503e+17,\n",
      "     -1.5282410383799979604e+16,\n",
      "     5.773449575496993589e+39,\n",
      "     5.3433937458969078745e+39,\n",
      "     469305,\n",
      "     0,\n",
      "     389174377867415552),\n",
      " 2: (2.1203328628257383575e+17,\n",
      "     64756017185446.597675,\n",
      "     1.7081708964304293988e+39,\n",
      "     1.3053204002923686526e+39,\n",
      "     222579,\n",
      "     56629159325661976,\n",
      "     317096453837818648),\n",
      " 3: (2.07697267279316528e+17,\n",
      "     2.6367481275468926898e+16,\n",
      "     8.87068088650291795e+38,\n",
      "     6.4929291690119722136e+38,\n",
      "     232617,\n",
      "     113258318651323952,\n",
      "     281862280268830256),\n",
      " 4: (5.5969475326659699072e+17,\n",
      "     1.33446813443680867414e+17,\n",
      "     9.2122867143093983015e+38,\n",
      "     9.627956722873028641e+38,\n",
      "     122946,\n",
      "     422111132843500776,\n",
      "     719865965107815656),\n",
      " 5: (6.456860630541701375e+17,\n",
      "     5.6144108286815278803e+17,\n",
      "     6.4715415135215808676e+40,\n",
      "     3.987460898575807787e+40,\n",
      "     1668640,\n",
      "     393134769365348352,\n",
      "     930014099572076544),\n",
      " 6: (5.495206401424045495e+17,\n",
      "     5.8497339784969710225e+17,\n",
      "     2.7716791620459931558e+40,\n",
      "     1.4838093329617353645e+40,\n",
      "     1250682,\n",
      "     324392439879838488,\n",
      "     737011256830966552),\n",
      " 7: (5.3164255361245846875e+17,\n",
      "     6.669296563733484617e+17,\n",
      "     4.5626733404635158524e+40,\n",
      "     2.1079075321851103023e+40,\n",
      "     2008965,\n",
      "     293103626956281392,\n",
      "     683953065157629488),\n",
      " 8: (8.9310204567759497444e+17,\n",
      "     4.968108793442648871e+17,\n",
      "     1.3543890113957935462e+39,\n",
      "     1.1327987709943823676e+39,\n",
      "     118980,\n",
      "     723251540325008616,\n",
      "     1040069742423074024),\n",
      " 9: (9.975517053850623995e+17,\n",
      "     1.1090356952206949738e+18,\n",
      "     6.3766381711546202313e+38,\n",
      "     4.9175748637866167454e+38,\n",
      "     344157,\n",
      "     934416127407980544,\n",
      "     1062614886965399552),\n",
      " 10: (8.988889746543214081e+17,\n",
      "      1.1136802629381279622e+18,\n",
      "      8.4472033903792019157e+39,\n",
      "      7.9521884857715649864e+39,\n",
      "      900915,\n",
      "      740431215345503000,\n",
      "      1063491898722519832),\n",
      " 11: (8.6822745952314770606e+17,\n",
      "      1.1882069866621482475e+18,\n",
      "      2.587423673737934511e+40,\n",
      "      2.6973534766185294057e+40,\n",
      "      2199288,\n",
      "      687266562658809392,\n",
      "      1063599953263306288)}\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from IPython.display import set_matplotlib_formats\n",
    "set_matplotlib_formats('svg')\n",
    "from matplotlib import pyplot as plt\n",
    "plt.rcParams['svg.fonttype'] = 'none'\n",
    "\n",
    "from glob import glob\n",
    "import numpy as np\n",
    "from pprint import pprint\n",
    "import tarfile\n",
    "from tqdm import tqdm_notebook\n",
    "\n",
    "# Klotz, Jerome H. \"UPDATING SIMPLE LINEAR REGRESSION.\"\n",
    "# Statistica Sinica 5, no. 1 (1995): 399-403.\n",
    "# http://www.jstor.org/stable/24305577\n",
    "def online_simple_regression(accumulator, x, y):\n",
    "    Ax_, Ay_, Sxy, Sxx, n_, minx, maxx = accumulator or (0, 0, 0, 0, 0, None, None)\n",
    "\n",
    "    first = n_ == 0\n",
    "    n = n_ + x.size\n",
    "    rt_n, rt_n_ = np.sqrt((n, n_), dtype=np.float128)\n",
    "\n",
    "    Ax = (Ax_*n_ + x.sum(dtype=np.float128))/n\n",
    "    Ay = (Ay_*n_ + y.sum(dtype=np.float128))/n\n",
    "    \n",
    "    minx = x.min() if first else min(minx, x.min())\n",
    "    maxx = x.max() if first else max(maxx, x.max())\n",
    "    \n",
    "    X = Ax if first else (Ax_*rt_n_ + Ax*rt_n)/(rt_n_ + rt_n)\n",
    "    Y = Ay if first else (Ay_*rt_n_ + Ay*rt_n)/(rt_n_ + rt_n)\n",
    "\n",
    "    Sxx += np.sum((x - X)**2)\n",
    "    Sxy += np.sum((x - X)*(y - Y))\n",
    "\n",
    "    return Ax, Ay, Sxy, Sxx, n, minx, maxx\n",
    "\n",
    "def conv_px(s):\n",
    "    w, h = s.split(b'x')\n",
    "    return int(w)*int(h)\n",
    "\n",
    "conv_fti = [b'I', b'P', b'B0', b'B1'].index\n",
    "\n",
    "def collect(filename, queues):\n",
    "    px, log_target_q, byte_size, frame_type = np.loadtxt(\n",
    "        filename, dtype=np.int64, delimiter=',',\n",
    "        converters={1: conv_px, 4: conv_fti},\n",
    "        skiprows=1, usecols=range(1, 5), unpack=True)\n",
    "\n",
    "    blog64q57_ibpp = np.round((\n",
    "        np.log2(px, dtype=np.float128) - np.log2(byte_size*8, dtype=np.float128)\n",
    "    )*2**57).astype(np.int64)\n",
    "    \n",
    "    # These are the fixed point found by repeating this whole process\n",
    "    boundaries = [\n",
    "        [0, 381625*2**40, 655352*2**40, 967797*2**40],\n",
    "        [0, 356802*2**40, 848173*2**40, 967797*2**40],\n",
    "        [0, 288436*2**40, 671307*2**40, 967797*2**40],\n",
    "        [0, 264708*2**40, 622760*2**40, 967797*2**40]\n",
    "    ]\n",
    "\n",
    "    for fti in np.unique(frame_type):\n",
    "        buckets = list(zip(boundaries[fti][:-1], boundaries[fti][1:]))\n",
    "        for bi, bucket in enumerate(buckets):\n",
    "            low, high = bucket\n",
    "            idx = (frame_type==fti) & (log_target_q >= low) & (log_target_q < high)\n",
    "            if np.sum(idx, dtype=int) == 0: continue\n",
    "            b = (bi << 2) | fti\n",
    "            x, y = log_target_q[idx], blog64q57_ibpp[idx]\n",
    "            queue = queues.get(b, ([], []))\n",
    "            queue[0].append(x)\n",
    "            queue[1].append(y)\n",
    "            queues[b] = queue\n",
    "\n",
    "def aggregate(queues, partials):\n",
    "    for b, queue in queues.items():\n",
    "        x, y = np.concatenate(queue[0]), np.concatenate(queue[1])\n",
    "        partials[b] = online_simple_regression(partials.get(b, None), x, y)\n",
    "    queues.clear()\n",
    "\n",
    "partials = dict()\n",
    "# https://ba.rr-dav.id.au/data/rav1e/rc-data.tar.xz\n",
    "with tarfile.open('rc-data.tar.xz', 'r:xz') as tf:\n",
    "    queues, last_name = dict(), None\n",
    "    for ti in tqdm_notebook(tf, total=1077*255, leave=False):\n",
    "        name = ti.name.split('/')[0]\n",
    "        if last_name and name != last_name:\n",
    "            aggregate(queues, partials)\n",
    "        last_name = name\n",
    "        collect(tf.extractfile(ti), queues)\n",
    "    aggregate(queues, partials)\n",
    "pprint(partials)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fixed-point approximation\n",
    "\n",
    "The regression results are converted to a fixed-point representation,\n",
    "with the exponent in Q6 and the scale in Q3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " I: exp=48 scale=36 bucket=0\n",
      " I: exp=61 scale=55 bucket=1\n",
      " I: exp=77 scale=129 bucket=2\n",
      " P: exp=69 scale=32 bucket=0\n",
      "B0: exp=84 scale=30 bucket=0\n",
      "B1: exp=87 scale=27 bucket=0\n",
      "B1: exp=139 scale=84 bucket=1\n",
      "B0: exp=120 scale=68 bucket=1\n",
      " P: exp=104 scale=84 bucket=1\n",
      "B1: exp=61 scale=1 bucket=2\n",
      "B0: exp=68 scale=4 bucket=2\n",
      " P: exp=83 scale=19 bucket=2\n"
     ]
    },
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg  PUBLIC '-//W3C//DTD SVG 1.1//EN'  'http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd'>\n",
       "<svg width=\"428pt\" height=\"361pt\" version=\"1.1\" viewBox=\"0 0 428 361\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<defs>\n",
       "<style type=\"text/css\">*{stroke-linecap:butt;stroke-linejoin:round;}</style>\n",
       "</defs>\n",
       "<path d=\"m0 361h428v-361h-428z\" fill=\"none\"/>\n",
       "<path d=\"m26.9 337h391v-326h-391z\" fill=\"#fff\"/>\n",
       "<path d=\"m128 337v-326\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-width=\".8\"/>\n",
       "<defs>\n",
       "<path id=\"e\" d=\"m0 0v3.5\" stroke=\"#000\" stroke-width=\".8\"/>\n",
       "</defs>\n",
       "<use x=\"128.481999\" y=\"336.86\" stroke=\"#000000\" stroke-width=\".8\" xlink:href=\"#e\"/>\n",
       "<text transform=\"rotate(-0 128 351)\" x=\"128.481999\" y=\"351.458437\" font-family=\"DejaVu Sans\" font-size=\"10px\" text-anchor=\"middle\">0</text>\n",
       "<path d=\"m383 337v-326\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-width=\".8\"/>\n",
       "<use x=\"382.520938\" y=\"336.86\" stroke=\"#000000\" stroke-width=\".8\" xlink:href=\"#e\"/>\n",
       "<text transform=\"rotate(-0 383 351)\" x=\"382.520938\" y=\"351.458437\" font-family=\"DejaVu Sans\" font-size=\"10px\" text-anchor=\"middle\">10</text>\n",
       "<path d=\"m77.7 337v-326\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<defs>\n",
       "<path id=\"b\" d=\"m0 0v2\" stroke=\"#000\" stroke-width=\".6\"/>\n",
       "</defs>\n",
       "<use x=\"77.674211\" y=\"336.86\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#b\"/>\n",
       "<path d=\"m179 337v-326\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<use x=\"179.289787\" y=\"336.86\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#b\"/>\n",
       "<path d=\"m230 337v-326\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<use x=\"230.097575\" y=\"336.86\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#b\"/>\n",
       "<path d=\"m281 337v-326\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<use x=\"280.905363\" y=\"336.86\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#b\"/>\n",
       "<path d=\"m332 337v-326\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<use x=\"331.713151\" y=\"336.86\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#b\"/>\n",
       "<path d=\"m26.9 268h391\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-width=\".8\"/>\n",
       "<defs>\n",
       "<path id=\"d\" d=\"m0 0h-3.5\" stroke=\"#000\" stroke-width=\".8\"/>\n",
       "</defs>\n",
       "<use x=\"26.925\" y=\"267.801483\" stroke=\"#000000\" stroke-width=\".8\" xlink:href=\"#d\"/>\n",
       "<text transform=\"rotate(-0 19.9 272)\" x=\"19.925\" y=\"271.600701\" font-family=\"DejaVu Sans\" font-size=\"10px\" text-anchor=\"end\">0</text>\n",
       "<path d=\"m26.9 13.8h391\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-width=\".8\"/>\n",
       "<use x=\"26.925\" y=\"13.762543\" stroke=\"#000000\" stroke-width=\".8\" xlink:href=\"#d\"/>\n",
       "<text transform=\"rotate(-0 19.9 17.6)\" x=\"19.925\" y=\"17.561762\" font-family=\"DejaVu Sans\" font-size=\"10px\" text-anchor=\"end\">10</text>\n",
       "<path d=\"m26.9 319h391\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<defs>\n",
       "<path id=\"c\" d=\"m0 0h-2\" stroke=\"#000\" stroke-width=\".6\"/>\n",
       "</defs>\n",
       "<use x=\"26.925\" y=\"318.60927\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#c\"/>\n",
       "<path d=\"m26.9 217h391\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<use x=\"26.925\" y=\"216.993695\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#c\"/>\n",
       "<path d=\"m26.9 166h391\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<use x=\"26.925\" y=\"166.185907\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#c\"/>\n",
       "<path d=\"m26.9 115h391\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<use x=\"26.925\" y=\"115.378119\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#c\"/>\n",
       "<path d=\"m26.9 64.6h391\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<use x=\"26.925\" y=\"64.570331\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#c\"/>\n",
       "<path d=\"m130 322 72.3-54.2 0.461-0.448 52.5-50 0.597-1.12 55.8-67.2\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#1f77b4\" stroke-linecap=\"square\" stroke-width=\"1.5\"/>\n",
       "<path d=\"m128 318 68.6-74 0.698-3 94.6-154 0.776-0.649 22.6-29.3\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#ff7f0e\" stroke-linecap=\"square\" stroke-width=\"1.5\"/>\n",
       "<path d=\"m138 304 45.9-60.3 1.29-4.39 72.7-136 0.603-1.55 56.9-60.5\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#2ca02c\" stroke-linecap=\"square\" stroke-width=\"1.5\"/>\n",
       "<path d=\"m148 286 29.7-40.4 1.98-3.82 68.9-150 0.584-3.17 66.3-63.2\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#d62728\" stroke-linecap=\"square\" stroke-width=\"1.5\"/>\n",
       "<path d=\"m26.9 337v-326\" fill=\"none\" stroke=\"#000\" stroke-linecap=\"square\" stroke-width=\".8\"/>\n",
       "<path d=\"m418 337v-326\" fill=\"none\" stroke=\"#000\" stroke-linecap=\"square\" stroke-width=\".8\"/>\n",
       "<path d=\"m26.9 337h391\" fill=\"none\" stroke=\"#000\" stroke-linecap=\"square\" stroke-width=\".8\"/>\n",
       "<path d=\"m26.9 10.7h391\" fill=\"none\" stroke=\"#000\" stroke-linecap=\"square\" stroke-width=\".8\"/>\n",
       "<path d=\"m365 77.4h45.2q2 0 2-2v-57.7q0-2-2-2h-45.2q-2 0-2 2v57.7q0 2 2 2z\" fill=\"#fff\" opacity=\".8\" stroke=\"#ccc\"/>\n",
       "<path d=\"m367 23.8h20\" fill=\"none\" stroke=\"#1f77b4\" stroke-linecap=\"square\" stroke-width=\"1.5\"/>\n",
       "<text transform=\"rotate(-0 395 27.3)\" x=\"395.301563\" y=\"27.298437\" font-family=\"DejaVu Sans\" font-size=\"10px\">I</text>\n",
       "<path d=\"m367 38.5h20\" fill=\"none\" stroke=\"#ff7f0e\" stroke-linecap=\"square\" stroke-width=\"1.5\"/>\n",
       "<text transform=\"rotate(-0 395 42)\" x=\"395.301563\" y=\"41.976562\" font-family=\"DejaVu Sans\" font-size=\"10px\">P</text>\n",
       "<path d=\"m367 53.2h20\" fill=\"none\" stroke=\"#2ca02c\" stroke-linecap=\"square\" stroke-width=\"1.5\"/>\n",
       "<text transform=\"rotate(-0 395 56.7)\" x=\"395.301563\" y=\"56.654687\" font-family=\"DejaVu Sans\" font-size=\"10px\">B0</text>\n",
       "<path d=\"m367 67.8h20\" fill=\"none\" stroke=\"#d62728\" stroke-linecap=\"square\" stroke-width=\"1.5\"/>\n",
       "<text transform=\"rotate(-0 395 71.3)\" x=\"395.301563\" y=\"71.332812\" font-family=\"DejaVu Sans\" font-size=\"10px\">B1</text>\n",
       "<defs>\n",
       "<clipPath id=\"a\">\n",
       "<rect x=\"26.9\" y=\"10.7\" width=\"391\" height=\"326\"/>\n",
       "</clipPath>\n",
       "</defs>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<Figure size 504x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(7, 6))\n",
    "plt.axis('equal')\n",
    "plt.xticks([0, 10])\n",
    "plt.yticks([0, 10])\n",
    "plt.minorticks_on()\n",
    "plt.grid(b=True, which='major')\n",
    "plt.grid(b=True, which='minor', alpha=0.2)\n",
    "\n",
    "segments = dict()\n",
    "for b, accumulator in partials.items():\n",
    "    Ax, Ay, Sxy, Sxx, n, minx, maxx = accumulator\n",
    "\n",
    "\n",
    "    fti = b & 3\n",
    "    beta = Sxy/Sxx\n",
    "    alpha = Ay - beta*Ax\n",
    "    exp = int(np.round(beta*2**6))\n",
    "    beta_ = exp/2**6\n",
    "    alpha_ = Ay - beta_*Ax\n",
    "    scale = int(np.round(np.exp2(3 - alpha_/2**57)))\n",
    "    label = ['I', 'P', 'B0', 'B1'][fti]\n",
    "    print('%2s: exp=%d scale=%d bucket=%d' % (label, exp, scale, b>>2))\n",
    "\n",
    "    xs, ys = segments.get(label, ([], []))\n",
    "    xs = [minx/2**57, maxx/2**57]\n",
    "    ys = [xs[0]*beta_ + alpha_/2**57, xs[1]*beta_ + alpha_/2**57]\n",
    "    xs_, ys_ = segments.get(label, ([], []))\n",
    "    xs_.extend(xs)\n",
    "    ys_.extend(ys)\n",
    "    segments[label] = (xs_, ys_)\n",
    "\n",
    "best = dict()\n",
    "for label, xy in segments.items():\n",
    "    plt.plot(xy[0], xy[1], label=label)\n",
    "    \n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The endpoints of each linear regression, rounding only the exponent, are detailed in the following output.\n",
    "We use a cubic interpolation of these points to adjust the segment boundaries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'B0': ([0.39294372842822706,\n",
      "         2.2002986504858395,\n",
      "         2.250924723555803,\n",
      "         5.114042917135398,\n",
      "         5.137773646423531,\n",
      "         7.379457452900416],\n",
      "        [-1.4148623981702452337,\n",
      "         0.95729093703037120824,\n",
      "         1.1300545385380504622,\n",
      "         6.4984011514997921106,\n",
      "         6.5594675085524934565,\n",
      "         8.941256552934184157]),\n",
      " 'B1': ([0.7858874568564541,\n",
      "         1.955812458298777,\n",
      "         2.0338149703000394,\n",
      "         4.745877754380938,\n",
      "         4.768869761992487,\n",
      "         7.380207231894769],\n",
      "        [-0.7078397903627942796,\n",
      "         0.88252700847286351155,\n",
      "         1.032874654152443289,\n",
      "         6.9231360133281455494,\n",
      "         7.0480332114869880963,\n",
      "         9.536964237487599543]),\n",
      " 'I': ([0.06648789020906937,\n",
      "        2.910837019748877,\n",
      "        2.928984366459142,\n",
      "        4.995073556916916,\n",
      "        5.018565704152715,\n",
      "        7.21693359533783],\n",
      "       [-2.1348326785650876126,\n",
      "        -0.0015708314102318150611,\n",
      "        0.016045673317236682116,\n",
      "        1.9852869329723028162,\n",
      "        2.0293453729387404098,\n",
      "        4.6742567420208319777]),\n",
      " 'P': ([0.0,\n",
      "        2.700439718141098,\n",
      "        2.72792045456319,\n",
      "        6.453269166068452,\n",
      "        6.483814370183836,\n",
      "        7.3733719613652795],\n",
      "       [-1.9815443653140106649,\n",
      "        0.92986720580686056843,\n",
      "        1.04808682998986787,\n",
      "        7.1017784861859188554,\n",
      "        7.1273316591101902764,\n",
      "        8.280976660173624982])}\n"
     ]
    }
   ],
   "source": [
    "pprint(segments)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Piecewise-linear fit\n",
    "\n",
    "We applied a 3-segment piecewise-linear fit. The boundaries were aligned to integer values of pixels-per-bit,\n",
    "while optimizing for similarity to a cubic interpolation of the control points\n",
    "(log-quantizer as a function of log-bitrate)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "I [1. 4.] [381625. 655352.]\n",
      "P [  2. 139.] [356802. 848173.]\n",
      "B0 [ 2. 92.] [288436. 671307.]\n",
      "B1 [  2. 126.] [264708. 622760.]\n"
     ]
    },
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg  PUBLIC '-//W3C//DTD SVG 1.1//EN'  'http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd'>\n",
       "<svg width=\"428pt\" height=\"361pt\" version=\"1.1\" viewBox=\"0 0 428 361\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<defs>\n",
       "<style type=\"text/css\">*{stroke-linecap:butt;stroke-linejoin:round;}</style>\n",
       "</defs>\n",
       "<path d=\"m0 361h428v-361h-428z\" fill=\"none\"/>\n",
       "<path d=\"m26.9 337h391v-326h-391z\" fill=\"#fff\"/>\n",
       "<path d=\"m128 337v-326\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-width=\".8\"/>\n",
       "<defs>\n",
       "<path id=\"f\" d=\"m0 0v3.5\" stroke=\"#000\" stroke-width=\".8\"/>\n",
       "</defs>\n",
       "<use x=\"128.482565\" y=\"336.86\" stroke=\"#000000\" stroke-width=\".8\" xlink:href=\"#f\"/>\n",
       "<text transform=\"rotate(-0 128 351)\" x=\"128.482565\" y=\"351.458437\" font-family=\"DejaVu Sans\" font-size=\"10px\" text-anchor=\"middle\">0</text>\n",
       "<path d=\"m383 337v-326\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-width=\".8\"/>\n",
       "<use x=\"382.521504\" y=\"336.86\" stroke=\"#000000\" stroke-width=\".8\" xlink:href=\"#f\"/>\n",
       "<text transform=\"rotate(-0 383 351)\" x=\"382.521504\" y=\"351.458437\" font-family=\"DejaVu Sans\" font-size=\"10px\" text-anchor=\"middle\">10</text>\n",
       "<path d=\"m77.7 337v-326\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<defs>\n",
       "<path id=\"d\" d=\"m0 0v2\" stroke=\"#000\" stroke-width=\".6\"/>\n",
       "</defs>\n",
       "<use x=\"77.674777\" y=\"336.86\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#d\"/>\n",
       "<path d=\"m179 337v-326\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<use x=\"179.290352\" y=\"336.86\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#d\"/>\n",
       "<path d=\"m230 337v-326\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<use x=\"230.09814\" y=\"336.86\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#d\"/>\n",
       "<path d=\"m281 337v-326\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<use x=\"280.905928\" y=\"336.86\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#d\"/>\n",
       "<path d=\"m332 337v-326\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<use x=\"331.713716\" y=\"336.86\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#d\"/>\n",
       "<path d=\"m26.9 268h391\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-width=\".8\"/>\n",
       "<defs>\n",
       "<path id=\"e\" d=\"m0 0h-3.5\" stroke=\"#000\" stroke-width=\".8\"/>\n",
       "</defs>\n",
       "<use x=\"26.925\" y=\"267.801483\" stroke=\"#000000\" stroke-width=\".8\" xlink:href=\"#e\"/>\n",
       "<text transform=\"rotate(-0 19.9 272)\" x=\"19.925\" y=\"271.600701\" font-family=\"DejaVu Sans\" font-size=\"10px\" text-anchor=\"end\">0</text>\n",
       "<path d=\"m26.9 13.8h391\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-width=\".8\"/>\n",
       "<use x=\"26.925\" y=\"13.762543\" stroke=\"#000000\" stroke-width=\".8\" xlink:href=\"#e\"/>\n",
       "<text transform=\"rotate(-0 19.9 17.6)\" x=\"19.925\" y=\"17.561762\" font-family=\"DejaVu Sans\" font-size=\"10px\" text-anchor=\"end\">10</text>\n",
       "<path d=\"m26.9 319h391\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<defs>\n",
       "<path id=\"c\" d=\"m0 0h-2\" stroke=\"#000\" stroke-width=\".6\"/>\n",
       "</defs>\n",
       "<use x=\"26.925\" y=\"318.60927\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#c\"/>\n",
       "<path d=\"m26.9 217h391\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<use x=\"26.925\" y=\"216.993695\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#c\"/>\n",
       "<path d=\"m26.9 166h391\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<use x=\"26.925\" y=\"166.185907\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#c\"/>\n",
       "<path d=\"m26.9 115h391\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<use x=\"26.925\" y=\"115.378119\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#c\"/>\n",
       "<path d=\"m26.9 64.6h391\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#b0b0b0\" stroke-linecap=\"square\" stroke-opacity=\".2\" stroke-width=\".8\"/>\n",
       "<use x=\"26.925\" y=\"64.570331\" stroke=\"#000000\" stroke-width=\".6\" xlink:href=\"#c\"/>\n",
       "<path d=\"m130 322 72.3-54.2 53.1-50.8 56.3-67.9\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#808080\" stroke-dasharray=\"3.7,1.6\"/>\n",
       "<defs>\n",
       "<path id=\"b\" d=\"m0 1.5c0.398 0 0.779-0.158 1.06-0.439 0.281-0.281 0.439-0.663 0.439-1.06s-0.158-0.779-0.439-1.06c-0.281-0.281-0.663-0.439-1.06-0.439s-0.779 0.158-1.06 0.439c-0.281 0.281-0.439 0.663-0.439 1.06s0.158 0.779 0.439 1.06c0.281 0.281 0.663 0.439 1.06 0.439z\" stroke=\"#808080\"/>\n",
       "</defs>\n",
       "<g clip-path=\"url(#a)\" fill=\"#808080\" stroke=\"#808080\">\n",
       "<use x=\"130.17083\" y=\"322.034545\" xlink:href=\"#b\"/>\n",
       "<use x=\"202.447779\" y=\"267.801483\" xlink:href=\"#b\"/>\n",
       "<use x=\"255.500521\" y=\"216.993695\" xlink:href=\"#b\"/>\n",
       "<use x=\"311.819845\" y=\"149.05716\" xlink:href=\"#b\"/>\n",
       "</g>\n",
       "<path d=\"m130 322 13.5-9.1 12.9-9.1 12.4-9.1 11.8-9.1 11.3-9.1 10.9-9.1 10.4-9.1 9.97-9.1 9.57-9.1 9.19-9.1 8.84-9.1 8.51-9.1 8.21-9.1 7.93-9.1 7.67-9.1 7.44-9.1 7.23-9.1 7.05-9.1 6.89-9.1\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#1f77b4\" stroke-linecap=\"square\" stroke-width=\"1.5\"/>\n",
       "<g stroke=\"#808080\">\n",
       "<path d=\"m128 318 69.1-75.7 95.2-155 22.9-29.5\" clip-path=\"url(#a)\" fill=\"none\" stroke-dasharray=\"3.7,1.6\"/>\n",
       "<g clip-path=\"url(#a)\" fill=\"#808080\">\n",
       "<use x=\"128.49652\" y=\"318.140425\" xlink:href=\"#b\"/>\n",
       "<use x=\"197.636603\" y=\"242.397589\" xlink:href=\"#b\"/>\n",
       "<use x=\"292.872367\" y=\"86.952659\" xlink:href=\"#b\"/>\n",
       "<use x=\"315.815404\" y=\"57.43243\" xlink:href=\"#b\"/>\n",
       "</g>\n",
       "</g>\n",
       "<path d=\"m128 318 15.2-13.7 13.9-13.7 12.7-13.7 11.6-13.7 10.6-13.7 9.84-13.7 9.16-13.7 8.61-13.7 8.2-13.7 7.91-13.7 7.75-13.7 7.72-13.7 7.82-13.7 8.06-13.7 8.42-13.7 8.91-13.7 9.53-13.7 10.3-13.7 11.2-13.7\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#ff7f0e\" stroke-linecap=\"square\" stroke-width=\"1.5\"/>\n",
       "<g stroke=\"#808080\">\n",
       "<path d=\"m138 304 45.9-61.3 74.2-140 57.4-61.4\" clip-path=\"url(#a)\" fill=\"none\" stroke-dasharray=\"3.7,1.6\"/>\n",
       "<g clip-path=\"url(#a)\" fill=\"#808080\">\n",
       "<use x=\"138.486838\" y=\"303.744497\" xlink:href=\"#b\"/>\n",
       "<use x=\"184.386199\" y=\"242.397589\" xlink:href=\"#b\"/>\n",
       "<use x=\"258.592894\" y=\"102.077606\" xlink:href=\"#b\"/>\n",
       "<use x=\"315.95348\" y=\"40.658749\" xlink:href=\"#b\"/>\n",
       "</g>\n",
       "</g>\n",
       "<path d=\"m138 304 12.5-13.8 11.1-13.8 9.88-13.8 8.84-13.8 7.99-13.8 7.32-13.8 6.84-13.8 6.55-13.8 6.45-13.8 6.53-13.8 6.8-13.8 7.26-13.8 7.9-13.8 8.74-13.8 9.76-13.8 11-13.8 12.4-13.8 13.9-13.8 15.7-13.8\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#2ca02c\" stroke-linecap=\"square\" stroke-width=\"1.5\"/>\n",
       "<g stroke=\"#808080\">\n",
       "<path d=\"m148 286 31.3-43.4 69.4-152 66.8-65\" clip-path=\"url(#a)\" fill=\"none\" stroke-dasharray=\"3.7,1.6\"/>\n",
       "<g clip-path=\"url(#a)\" fill=\"#808080\">\n",
       "<use x=\"148.440733\" y=\"285.783369\" xlink:href=\"#b\"/>\n",
       "<use x=\"179.787342\" y=\"242.397589\" xlink:href=\"#b\"/>\n",
       "<use x=\"249.183659\" y=\"90.551403\" xlink:href=\"#b\"/>\n",
       "<use x=\"315.947811\" y=\"25.525455\" xlink:href=\"#b\"/>\n",
       "</g>\n",
       "</g>\n",
       "<path d=\"m148 286 11.6-13.7 9.93-13.7 8.53-13.7 7.35-13.7 6.41-13.7 5.71-13.7 5.25-13.7 5.02-13.7 5.03-13.7 5.27-13.7 5.75-13.7 6.47-13.7 7.42-13.7 8.61-13.7 10-13.7 11.7-13.7 13.6-13.7 15.7-13.7 18.1-13.7\" clip-path=\"url(#a)\" fill=\"none\" stroke=\"#d62728\" stroke-linecap=\"square\" stroke-width=\"1.5\"/>\n",
       "<path d=\"m26.9 337v-326\" fill=\"none\" stroke=\"#000\" stroke-linecap=\"square\" stroke-width=\".8\"/>\n",
       "<path d=\"m418 337v-326\" fill=\"none\" stroke=\"#000\" stroke-linecap=\"square\" stroke-width=\".8\"/>\n",
       "<path d=\"m26.9 337h391\" fill=\"none\" stroke=\"#000\" stroke-linecap=\"square\" stroke-width=\".8\"/>\n",
       "<path d=\"m26.9 10.7h391\" fill=\"none\" stroke=\"#000\" stroke-linecap=\"square\" stroke-width=\".8\"/>\n",
       "<path d=\"m365 77.4h45.2q2 0 2-2v-57.7q0-2-2-2h-45.2q-2 0-2 2v57.7q0 2 2 2z\" fill=\"#fff\" opacity=\".8\" stroke=\"#ccc\"/>\n",
       "<path d=\"m367 23.8h20\" fill=\"none\" stroke=\"#1f77b4\" stroke-linecap=\"square\" stroke-width=\"1.5\"/>\n",
       "<text transform=\"rotate(-0 395 27.3)\" x=\"395.301563\" y=\"27.298437\" font-family=\"DejaVu Sans\" font-size=\"10px\">I</text>\n",
       "<path d=\"m367 38.5h20\" fill=\"none\" stroke=\"#ff7f0e\" stroke-linecap=\"square\" stroke-width=\"1.5\"/>\n",
       "<text transform=\"rotate(-0 395 42)\" x=\"395.301563\" y=\"41.976562\" font-family=\"DejaVu Sans\" font-size=\"10px\">P</text>\n",
       "<path d=\"m367 53.2h20\" fill=\"none\" stroke=\"#2ca02c\" stroke-linecap=\"square\" stroke-width=\"1.5\"/>\n",
       "<text transform=\"rotate(-0 395 56.7)\" x=\"395.301563\" y=\"56.654687\" font-family=\"DejaVu Sans\" font-size=\"10px\">B0</text>\n",
       "<path d=\"m367 67.8h20\" fill=\"none\" stroke=\"#d62728\" stroke-linecap=\"square\" stroke-width=\"1.5\"/>\n",
       "<text transform=\"rotate(-0 395 71.3)\" x=\"395.301563\" y=\"71.332812\" font-family=\"DejaVu Sans\" font-size=\"10px\">B1</text>\n",
       "<defs>\n",
       "<clipPath id=\"a\">\n",
       "<rect x=\"26.9\" y=\"10.7\" width=\"391\" height=\"326\"/>\n",
       "</clipPath>\n",
       "</defs>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<Figure size 504x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(7, 6))\n",
    "plt.axis('equal')\n",
    "plt.xticks([0, 10])\n",
    "plt.yticks([0, 10])\n",
    "plt.minorticks_on()\n",
    "plt.grid(b=True, which='major')\n",
    "plt.grid(b=True, which='minor', alpha=0.2)\n",
    "\n",
    "from scipy import optimize\n",
    "\n",
    "for ft, xy in segments.items():\n",
    "    f = np.poly1d(np.polyfit(np.array(xy[1]).astype(float), np.array(xy[0]).astype(float), 3))\n",
    "    ys = np.linspace(min(xy[1]), max(xy[1]), 20)\n",
    "    def cost(X):\n",
    "        y0 = np.array([ys[0], X[0], X[1], ys[-1]]).astype(float)\n",
    "        x0 = f(y0)\n",
    "        f0 = np.where(ys<X[0],\n",
    "                      np.poly1d(np.polyfit(y0[:2], x0[:2], 1))(ys),\n",
    "                      np.where(ys<X[1],\n",
    "                              np.poly1d(np.polyfit(y0[1:3], x0[1:3], 1))(ys),\n",
    "                              np.poly1d(np.polyfit(y0[2:], x0[2:], 1))(ys)))\n",
    "        return ((f0-f(ys))**2).sum()\n",
    "    X = optimize.fmin(cost, [2, 5], disp=0)\n",
    "    X = np.log2(np.ceil(np.exp2(X)))\n",
    "    print(ft, np.exp2(X), np.round(f(X)*2**17))\n",
    "    y0 = [ys.min(), X[0], X[1], ys.max()]\n",
    "    x0 = f(y0)\n",
    "    plt.plot(x0, y0, '.--', lw=1, c='grey')\n",
    "    plt.plot(f(ys), ys, label=ft)\n",
    "\n",
    "plt.legend();"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
